% run this file to recreate Figure 8 ("The Impact of Oil Inventores")

% The figure includes the IRFs for 4 different cases: i) with and without oil storage for oil-supply shocks
% and ii) with and without oil storage for TFP shocks

comp = [0 1];
track = 0;

% IRFs for full information
% Calibration of oil supply shocks, early 2000s
rhoe_p = 0.99999;   % serial correlation of persistent energy shocks
sige_p =  0.0655*0.01; % std. deviation of innovations to those shocks
rhoe_t = 0.8158;   % serial correlation of transitory energy shocks
sige_t = 0.01;    % std. deviation of innovations to those shocks

%Calibration for TFP shocks "early 2000s"
rhoz_p = 0.99999;   % serial correlation of persistent productivity shocks
sigz_p = 0.2707*0.3694*0.01; % std. deviation of innovations to those shocks
rhoz_t = 0.8466;   % serial correlation of transitory productivity shocks
sigz_t = 0.3694*0.01;    % std. deviation of innovations to those shocks

whatinfo = 2;
% whatinfo = 2 implies "learning" solution 
answinv = 1;
% answinv = 1 implies solution where oil inventoires are present  
answlearn = 1;
% answlearn = 1 indicates the usual learning solution. "answlearn" = 0 implies the learning solution used to construct Table 5 

param = [rhoe_p sige_p rhoe_t sige_t rhoz_p sigz_p rhoz_t sigz_t answinv whatinfo answlearn];

gainkeep = []; %this is not used here but in other programs





%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 1- Compute the IRFs when oil inventories are present

% Call the program that computes the s.s. and the linearization

lin
system = 'lin';

% Call the KW(2002) program

[M,PI,G]=resolkw(system,track,comp,param,gainkeep);

% compute 40-periods IRFs for all variables following a permanent oil-supply shock
nM = rows(M);
nPI = rows(PI);
S = zeros(1,nM)';
S(6) = -0.01271; % Shock magnitude such that the full-information impact response of oil prices is 10%
nir = 50;
IRB = zeros(nir,nPI+1);




for i=1:nir
   %Index of Period;
   IRB(i,1) = i-1;
   %Observables
   IRB(i,2:nPI+1)=(PI*S)';
   % Update the state vector:;
   S = M*S;  
end
IR1 = IRB; % save responses

%now compute responses for shock to TFP
S = zeros(1,nM)';
S(8) = 0.0390; % Shock magnitude such that the full-information impact response of oil prices is 10%
nir = 50;
IRB = zeros(nir,nPI+1);

for i=1:nir
   %Index of Period;
   IRB(i,1) = i-1;
   %Observables
   IRB(i,2:nPI+1)=(PI*S)';
   % Update the state vector:;
   S = M*S;  
end
IR2 = IRB; % save responses
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 2- Compute IRS for when oil inventoires are not present.

whatinfo = 2;
answinv = 0;
% answinv = 1 implies solution where oil inventoires are present 

param = [rhoe_p sige_p rhoe_t sige_t rhoz_p sigz_p rhoz_t sigz_t answinv whatinfo answlearn];


lin
system = 'lin';

% Call the KW(2002) program

[M,PI,G]=resolkw(system,track,comp,param,gainkeep);

% compute IRFs for all variables following a permanent oil-supply shock
nM = rows(M);
nPI = rows(PI);
S = zeros(1,nM)';
S(6) = -0.01240; % Shock magnitude such that the full-information impact response of oil prices is 10%
nir = 50;
IRB = zeros(nir,nPI+1);

for i=1:nir
   %Index of Period;
   IRB(i,1) = i-1;
   %Observables
   IRB(i,2:nPI+1)=(PI*S)';
   % Update the state vector:;
   S = M*S;  
end
IR3 = IRB; % save responses

%now compute responses for the shocks to TFP
S = zeros(1,nM)';
S(8) = 0.0442; % Shock magnitude such that the full-information impact response of oil prices is 10%
nir = 50;
IRB = zeros(nir,nPI+1);


for i=1:nir
   %Index of Period;
   IRB(i,1) = i-1;
   %Observables
   IRB(i,2:nPI+1)=(PI*S)';
   % Update the state vector:;
   S = M*S;  
end
IR4 = IRB; % save responses 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%




%Express all IRF is percentage points for graphiing purposes
IR1 = 100*IR1;
IR2 = 100*IR2;
IR3 = 100*IR3;
IR4 = 100*IR4;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



% Now ready to print responses

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(8)
% a 3 by 2 graph with oil-supply (left panel) and TPF (right panel) shocks side
% by side





myylabel = 'Deviation from s.s. (%)  ';
mysize = 15;
myxlabel = 'Quarters';

fmt = 5.3;
epais = 4.5;


TT2 = 12;
dt2 = (0:TT2-1)';

trait1 = 'r-';
trait2 = 'g:';


subplot(321),plot(dt2,IR1(1:TT2,1+iy),trait1,dt2,IR3(1:TT2,1+iy),trait2,'LineWidth',epais),
title('GDP','FontWeight','bold','FontSize',mysize)
%axis([0 TT2-1 -0.4 -0.1])
ylabel(myylabel,'FontWeight','bold','FontSize',mysize)
set(gca,'FontSize',mysize);
xt = get(gca, 'YTick');
xt1 = fortick(xt,fmt);
set(gca, 'YTickLabel', xt1, 'YTick', xt);


h1 = text(3.0,-0.04,'Oil Supply Shocks');
h1.Interpreter = 'latex';
h1.FontWeight = 'bold';
h1.FontAngle = 'italic';
h1.FontSize = 20;



subplot(323),plot(dt2,IR1(1:TT2,1+is),trait1,dt2,IR3(1:TT2,1+is),trait2,'LineWidth',epais),
%axis([0 TT2-1 -8.0 0.0])
title('Oil Inventories','FontWeight','bold','FontSize',mysize)
ylabel(myylabel,'FontWeight','bold','FontSize',mysize)
set(gca,'FontSize',mysize);
xt = get(gca, 'YTick');
xt1 = fortick(xt,fmt);
set(gca, 'YTickLabel', xt1, 'YTick', xt);


subplot(325),plot(dt2,IR1(1:TT2,1+ipo),trait1,dt2,IR3(1:TT2,1+ipo),trait2,'LineWidth',epais),
%axis([0 TT2-1  4.0 11.0])
title('Oil Prices','FontWeight','bold','FontSize',mysize)
xlabel(myxlabel,'FontSize',mysize,'FontWeight','bold')
ylabel(myylabel,'FontWeight','bold','FontSize',mysize)
set(gca,'FontSize',mysize);
xt = get(gca, 'YTick');
xt1 = fortick(xt,fmt);
set(gca, 'YTickLabel', xt1, 'YTick', xt);



subplot(322),plot(dt2,IR2(1:TT2,1+iy),trait1,dt2,IR4(1:TT2,1+iy),trait2,'LineWidth',epais),
title('GDP','FontWeight','bold','FontSize',mysize)
%axis([0 TT2-1 -0.4 -0.1])
ylabel(myylabel,'FontWeight','bold','FontSize',mysize)
set(gca,'FontSize',mysize);
xt = get(gca, 'YTick');
xt1 = fortick(xt,fmt);
set(gca, 'YTickLabel', xt1, 'YTick', xt);


h1 = text(3.5,5.2,'TFP Shocks');
h1.Interpreter = 'latex';
h1.FontWeight = 'bold';
h1.FontAngle = 'italic';
h1.FontSize = 20;



subplot(324),plot(dt2,IR2(1:TT2,1+is),trait1,dt2,IR4(1:TT2,1+is),trait2,'LineWidth',epais),
%axis([0 TT2-1 -8.0 0.0])
title('Oil Inventories','FontWeight','bold','FontSize',mysize)
ylabel(myylabel,'FontWeight','bold','FontSize',mysize)
set(gca,'FontSize',mysize);
xt = get(gca, 'YTick');
xt1 = fortick(xt,fmt);
set(gca, 'YTickLabel', xt1, 'YTick', xt);


subplot(326),plot(dt2,IR2(1:TT2,1+ipo),trait1,dt2,IR4(1:TT2,1+ipo),trait2,'LineWidth',epais),
%axis([0 TT2-1  4.0 11.0])
title('Oil Prices','FontWeight','bold','FontSize',mysize)
xlabel(myxlabel,'FontSize',mysize,'FontWeight','bold')
ylabel(myylabel,'FontWeight','bold','FontSize',mysize)
set(gca,'FontSize',mysize);
xt = get(gca, 'YTick');
xt1 = fortick(xt,fmt);
set(gca, 'YTickLabel', xt1, 'YTick', xt);









legende1 = 'With Oil Storage';
legende2 = 'Without Oil Storage';
h = legend(legende1,legende2);
h.Interpreter = 'latex';
h.FontWeight = 'bold';
h.FontSize = 20;
h.Orientation = 'horizontal';










